function[G] = numgrad(FUN,x0,varargin)

% PURPOSE: computes numerical Jacobian matrix
%
% INPUTS: 

%     FUN = Handle or name for function that generates a n x 1 vector of %	individual LLH components
%     Note: That function is specific to a given application & should be 
%         generated & saved w/in that application
%         x0 = k x 1 vector of betas
%         varargin = all inputs needed for evaluating FUN, other than x0    
%         
%         
% OUTPUT: G = n x k numerical Jacobian
%
%   Written by:
%   Klaus Moeltner
%   Department of Applied Economics & Statistics
%   University of Nevada, Reno
%   based on GAUSS' "gradp" command


L0 = feval(FUN,x0,varargin{:}); % n by 1
n=length(L0);
k = length(x0);

G=zeros(n,k);

% Compute step size (amount of perturbation of x0)

ax0 = abs(x0);
dax0 = x0./ax0;
dh=(1e-8)*max([ax0 1e-2*ones(k,1)]')'.*dax0; % k x 1
xdh = x0+dh;
dh = xdh-x0; % apparently increases precision

int1 = diag(dh); % k by k square matrix with xdh on diagonal
int2 = repmat(x0,1,k); % k by k
arg = int2+int1;

i=1;
for i=1:k;
    G(:,i)=feval(FUN,arg(:,i),varargin{:});
end

int1=repmat(L0,1,k); % n by k
int2=repmat(dh',n,1);
G = (G-int1)./int2;
